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We analyze a complex scalar field with (f> interaction and a chem- 
ical potential \x on the lattice. An exact flux representation of the 
partition sum is used which avoids the complex action problem and 
based on a generalized worm algorithm we can run Monte Carlo simu- 
lations at arbitrary densities. We study thermodynamical quantities as 
a function of the chemical potential fj, for zero- and finite temperature. 
It is shown that at zero temperature thermodynamical observables are 
independent of [i up to a critical value /x c (Silver Blaze phenomenon). 
In a spectroscopy calculation we cross-check that fx c agrees with the 
mass m of the scalar field. The Silver Blaze region ends in a second 
order phase transition and we show that for low temperatures the sec- 
ond order phase boundary persists and separates a pseudo Silver Blaze 
region from a condensed phase with strong /z-dependence. 
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1 Introduction 



The QCD phase diagram is a topic which currently sees a lot of attention 
in experimental and theoretical work. A particularly interesting (and rather 
unexplored) part is the region at low temperatures and medium to high 
densities (medium to high values of the chemical potential n). While for 
large chemical potential exotic states of matter are conjectured, for small 
chemical potential and vanishing temperature one expects physics to be 
independent of \i: As long as the chemical potential is smaller than the 
mass of the state it couples to no excitations can be generated and physics 
must be independent of \x. This is an interesting non-perturbative effect 
which has come to be known as Silver Blaze problem [T]. 

In principle lattice QCD would be an adequate method to study such 
non-perturbative phenomena. However, when coupling a chemical potential 
the fermion determinant becomes complex and cannot be used as a prob- 
ability weight in a Monte Carlo simulation. This so called complex action 
problem is a common feature in many field theories at non-zero \x and has 
considerably held up the study of QCD and other systems at finite den- 
sity. Using new ideas such as worm algorithms [2], designed for updating 
systems with local constraints, several interesting models at finite density 
were rewritten in terms of a flux representation where the complex action 
problem is avoided and worm-based Monte Carlo simulations are possible 
(see, e.g., HI El El OH] for examples). This opens the possibility to 
explore the Silver Blaze phenomenon in an ab-initio lattice study. 

In this article we analyze the Silver Blaze problem for a charged scalar 
4 field. For scalar fields a suitable flux representation [TJ [10] can be 
found which avoids the complex action problem at finite fi and Monte Carlo 
simulations at high densities become possible. The model has been analyzed 
previously using the complex Langevin approach and mean field techniques 
|1 11 IT2] , In this paper we systematically study the transition at the end 
of the Silver Blaze region at zero temperature and attempt to follow the 
transition line for non-zero temperatures. This work sheds light on Silver 
Blaze type of phenomena and explores new strategies for overcoming the 
complex action problem in lattice calculations at finite density. 

The article is organized as follows: In the next section we present the 
flux representation and our observables. In Section 3 the worm algorithm 
is described (see also the appendix for a pseudo-code listing) and tested. 
The numerical analysis at vanishing temperature is discussed in Section 4, 
followed by the results for non-zero temperature in Section 5. A summary 
and discussion section completes the paper. 
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2 Flux representation for the charged scalar field 



In the conventional representation the lattice action of the charged scalar 
field with a 4 -interaction is given by 




+ A|<fe| 4 - [e"*- 4 ^x+s + e-^tifa-v] J . (1) 



The first sum is over the sites x of a N% x N t lattice with periodic boundary 
conditions, the second sum is over the four Euclidean directions v = 1, 2, 3, 4, 
and denotes the unit vector in ^-direction. The degrees of freedom are the 
complex valued field variables (p x attached to the sites x of the lattice, rj 
denotes the combination 8 + m 2 , where m is the bare mass parameter. The 
(f) A coupling is denoted by A and the chemical potential by jjl. All parameters 
are in units of the lattice spacing a, in other words the lattice spacing is set 
to a = 1 throughout this paper, except for some places where we display 
a explicitly for clarity. The partition sum Z is obtained by integrating the 
Boltzmann factor e~ s over all field configurations, Z = f D[(f>]e~ s . The 
measure is a product over all lattice sites x, with (f> x being integrated over 
the complex plane, i.e., D[(f>] = Y\ x J € d4> x /2iT. The normalization with 2-7T 
will be useful later. 

Note that in its present form the model has a complex action problem: 
For ^ / the sum of the temporal nearest neighbor terms is complex. 
Thus the Boltzmann factor cannot be used as a probability weight in a 
Monte Carlo simulation. The problem can be overcome by using alternative 
degrees of freedom, i.e., a flux representation. The first step to obtain the 
flux representation of the partition sum is an expansion of the Boltzmann 
factors for the nearest neighbor terms 

nexp(e^ 4 ## x+ s) exp(e-^ 4 <^ + ~) = (2) 



1 \ TT .1 \ / TT/ , \ n 
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where by ^{nn} we denote the sum over all configurations of the expansion 
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variables n XtU ,n XtU G [0, oo), i.e., the multiple sum 



oo oo 



E - n E E <»> 

{n,n} x,v n x , v =0 n x , v =0 

The complex field variables are now written in polar form, <p x = r x e lSx . 
Rewriting also the integration measure in polar coordinates one finds for 
the partition sum 

z = ^2 ( n 1 HIT r d0x c~ i6xT - jnx - v ~^ x - v ~ (n ^- v -"~^- v -"A 

{n,n} \x,u nx ^^ l J\x J— 27F / 

x (IT e^'" 1 ' 4- " 1 ' 4 ' / dr x r x + ^^ nx ' v+nx ~ v '" +nx ' v+nx ~ v '^ e~ VVx ~ XVx \ . (4) 



The integrals over the phase give rise to Kronecker deltas, which for nota- 
tional convenience we will write as <5(n). The integrals over the modulus we 
abbreviate as 

/•oo 

W(n) = / drr^e'^ ~ Xr . (5) 
Jo 

They can easily be computed numerically. The partition sum now reads: 

z = E I L L , II s ( E - - ( n *-v - n x-v, v )] 

x ^JJ e »ln xA -n xA ] w ^ ^ + ^ v + ^ ^ + nx _ p ^ ^ (6) 

In this form the complex phase problem is gone. All weight factors for 
configurations of the n and n variables are real and non-negative. However, 
configurations of the n and n can still have weight 0. This is always the case 
if at a site x the constraint enforced by the Kronecker deltas is violated, i.e., 
when at least for one x one has ^2i V \n x ^ v — n Xjl/ — [n x -v,v — n x --p,v)] ^ 0. 

In the current representation the constraints mix both, the n and the n 
variables. The structure of the constraints can be simplified by introducing 
new variables fc Xj „ G (— oo, oo) and l XjU G [0, oo). They are related to the old 
variables by 

i^x,u ^x,v — k Xjl/ and ^x,v ~\~ f^x,v — ~\~ ^x,v • (7) 
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The partition sum turns into 



z = E n ,, ,/,,„.,, n« st^-* 




x—u,v 



+ 2(l XtV + l x -p, v )]jj. (8) 

In the final form (|8j), which we now refer to as flux representation, the 
constraints no longer mix the two types of flux variables. Obviously only 
the fc-fluxes are subject to conservation of flux at each site x, i.e., only they 
must obey ^2 u [k XjU — k x ^ v ] = for all x. 

Also observables need to be mapped to the flux representation. For 
thermodynamic observables which are obtained as derivatives of In Z with 
respect to the various couplings this is particularly simple. We will need the 
following observables which we refer to as particle number density n, the 
particle number susceptibility Xn, the derivative x' n °f Xn with respect to fx, 
the field expectation value (|</> 2 |), the corresponding susceptibility x\<f>\ 2 an d 
its derivative xUvz with respect to fi. They are defined as 

_ TdlnZ 1 dlnZ dn , _ d X n 

n ~ v~d~iT ~ N[W t dy. ' Xn ~ ' Xn ~ ~dJT ' 

-T dlnZ __ -1 dlnZ _ -d{\(j)\ 2 ) , dx\w 

V drj Nf N t dr, ' XW ~ drj ' X ^ 2 ~ dy ' 

All observables are normalized by the spatial volume V = Ng in order to 
make them intensive quantities. We remark, that the observables ^ are 
dimensionful (except X\<p\ 2 )- Since we here use a fixed scale approach (the 
lattice spacing a is set by r] and A and kept fixed) and drive the temperature 
T = 1/aNt by changing Nt, we will often show results in lattice units. Where 
a comparison of results at different scales is necessary we make quantities 
dimensionless using suitable powers of fj, c , the critical value of the chemical 
potential at zero temperature (see the discussion below). 

It is straightforward to work out the derivatives that define the observ- 
ables Q also for the flux representation of Z and we only discuss two ex- 
amples to illustrate the general structure. Particularly simple is the particle 
number density. One immediately sees that the flux representation of n is 
given by 
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i.e., the particle number is simply the net flux of k around compactified time 
(= the 4-direction) . The expectation value (...) is now of course understood 
as a vacuum expectation value in the flux representation. Not much more 
difficult is the flux representation for (\4>\ 2 ) = —djdr\ In Z/NfNt- Note that 
77 enters only through the weight factors W(n) defined in ^ which have the 
obvious property —d/dr]W(n) = W(n + 2). Application of the product rule 
to the product of the factors W in ^ and subsequent completion of the 
remaining products to the full Boltzmann factor gives rise to 



where we use the abbreviation f x = Yl v [\kx,u\ + \k x -„^\ + 2(l X]V + l x -tr >v )]. 
In a similar way the other observables listed in ([9]) can be obtained as higher 
moments of the fluxes and the sum of ratios of W. 

3 Worm algorithm and tests 

Having discussed the flux representation of the partition sum and the ob- 
servables, we now address the task of finding a suitable algorithm. In the 
final representation Q of the partition sum we use two different sets of 
variables. One of them, the variables l xu , are not subject to a constraint. 
For these we simply apply sweeps of local Metropolis updates: Each link 
(x, v) is visited and we generate a trial configuration by increasing or de- 
creasing l XjV by ±1 with equal probability. The change is accepted with the 
usual Metropolis probability given by the ratio of the weights of new and 
old configuration. A negative trial link is always rejected. 

Considerably more interesting is the update of the flux variables k X}V 
where a suitable update needs to take into account the local constraints 
Ylv[kx,v — k x _-p :U ] = 0. For that purpose we use a generalization of the 
Prokof'ev-Svistunov worm algorithm [2]. The worm starts at a randomly 
chosen site of the lattice and then grows by adding links to form a chain of 
links. At each position x the worm chooses randomly the next link from the 
8 links attached to x. For that link it proposes to change the corresponding 
variable k by ±1 and accepts that change with a local Metropolis step. At 
the starting point and the head of the worm the constraint is violated, and 
the worm must continue until its head reaches the starting point such that 
along the closed contour of changed links all constraints are intact again. 
A possible initial configuration for the k xv is to set k xv = at all links, a 
configuration that clearly obeys the constraints. 




(11) 
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The situation for a worm update of the variables k x ,u is a little bit tricky: 
The local weight for a variable k X)l/ not only has factors that live on the link 
(the factorials in i/| + lx,uY- and the contribution from the chemical 

potential), but also two factors at the two endpoints of the link, the terms 
W(f x ) and W(f x+ -p) (where again f x = YlA\ k x,u\MK-V,v\+2(lx,v+lx--D,v)])- 
The contributions from the factors W at the endpoints need to be suitably 
distributed among subsequent steps of the worm. In the appendix we provide 
pseudo-code for the initialization of the worm and the subsequent steps, 
illustrating how we solved the weight assignment for the W contributions. 
Worms for the variables k x ,u are alternated with the local Metropolis sweeps 
for the unconstrained variables l x>v . 

In our simulations we use Nf x Nt lattices with N s ranging from 4 to 
24, and Nt from 2 to 100. A full update sweep is the combination of one 
sweep for the l x v variables and one worm for the k x ,v fluxes. We typically 
use between 500.000 and 1.000.000 measurements separated by 5 full update 
sweeps. The number of equilibration sweeps is typically 25.000. All errors 
we quote are statistical errors determined with the jackknife method. The 
weig hts W defined in ^ were evaluated numerically using Mathematica and 
pre-stored for the Monte Carlo simulation. 

The algorithm was thoroughly tested by comparing it to the free case 
(A = 0) where exact results can be obtained using Fourier transformation, 
and to a conventional simulation at \x = 0. In Fig. [T] we show examples for 
the comparison using 8 4 lattices. The lhs. plot shows the particle number 
density n as a function of 77 for different values of the chemical potential fx 
for the case of A = 0. The symbols are the data from our flux simulation 
and the full curves the results from Fourier transformation. The rhs. plot 
compares the fj, = results for the expectation value (|</>| 2 ) from the flux 
simulation (triangles) to those of a conventional simulation (crosses). We 
plot (\4>\ 2 ) as a function of r\ for different values of A. The two plots show 
that the results from the flux simulation very nicely match the exact curves 
for the free case and the data from a conventional simulation at fj, = 0. 
Further comparisons were done for other observables and more values of the 
parameters - all with very good agreement. 

Finally we performed a comparison with the complex Langevin results 
at finite \i by Aarts [TT] , and we found excellent agreement of our data with 
the complex Langevin results that are available on lattices of sizes 4 4 to 
10 4 . It is interesting to note that these results are very well described also 
by mean- field theory [12* . In that study the quartic coupling was fixed to 
A = 1.0, and the boson mass M satisfied (aM) 2 = 1.47. As a result of 
this heavy mass, the temperatures reached on lattices L 4 = 4 4 to 10 4 are 
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Figure 1: Lhs. plot: Comparison of results for the particle number density n 
from the flux simulation at A = (symbols) to the exact result obtained from 
Fourier transformation (curves). We plot n as a function of ij for different 
values of /i using lattices of size 8 4 . Rhs. plot: Results for the expectation 
value (\4>\ 2 ) as a function of rj for different values of A, again on lattices of 
size 8 4 . The chemical potential is set to [i = and we compare the results 
from the flux simulation (triangles, circles, squares) to the outcome of a 
conventional [i = simulation (crosses). 

in the range T/M = 0.08, .., 0.21. At the same time, the system is large, 
since ML = 5, .., 12. So the observed mean- field behavior should come as 
no surprise. Our comparison with [TTj [T2] confirms the findings by Aarts. 

To summarize: Our tests establish that the flux representation is correct, 
that the algorithm works, and that the observables were correctly mapped 
to the corresponding expressions in terms of the fluxes. 

4 Transition at zero temperature 

We begin our study of the phase diagram with an analysis of the situation at 
zero temperature. For this analysis we set the temporal extent of the lattice 
to Nt = 100, such that in lattice units the temperature is T = 1/Nt = 0.01. 
Although not exactly 0, this is much smaller than all other involved scales. 
To study the approach of our observables to the thermodynamic limit we 
analyze different spatial volumes using N s = 4, 8, 12, 16, 20 and 24. We 
consider two different combinations of the couplings r\ and A. The two points 
in parameter space are r\ = 9.0, A = 1.0 and rj = 7.44, A = 1.0. 
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Figure 2: Results at T = 0.01 for ij = 9.0 and A = 1.0. In the lhs. column 
of plots we show n, \n and x' n as a function of /i (top to bottom). In the 
rhs. column we show (|</>| 2 ), X\4>\ 2 an d Xup- We compare the results for 
4 3 x 100, 8 3 x 100, 12 3 x 100, 16 3 x 100 and 20 3 x 100. 
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In Fig. [| we show the rj = 9.0, A = 1.0 results for n, Xn and x'n as a 
function of ji in the lhs. column of plots, and for (|0| 2 ), X\<f>\ 2 an( i Xu\a on 
the rhs. The particle number density n vanishes up to a critical value of 
He = 1.146(1) which marks the end of the Silver Blaze region where observ- 
ables are independent of ji. For values \i > ji c the particle number density 
rises essentially linearly. Consequently the particle number susceptibility Xn 
is discontinuous at /x c , jumping from Xn = to a finite value at ji c - The 
comparison of the different spatial volumes illustrates how the discontinuity 
emerges in the thermodynamic limit. Finally, when looking at the third 
/i-derivative xL °f the free energy we observe that for increasing volume a 
strongly peaked maximum emerges that very precisely marks the value of 
fj, c = 1.146(1). Furthermore, the height of the peak grows with increasing 
volume V = Ng (for a quantitative assessment see below), indicating the 
emergence of a true singularity in the thermodynamic limit. The observ- 
ables (|</>| 2 ), X\<f>\ 2 an d Xuia behave similarly to the derivatives with respect 
to ji: The expectation value (\4>\ 2 ) is constant below ji c and for \i > \i c 
increases with constant slope. The corresponding susceptibility X\M 2 thus 
also develops a step in the thermodynamic limit. Finally, another derivative 
with respect to ji, i.e., Xui2> shows the corresponding peaks, which again 
grow with the volume. 

The results for the second parameter set, rj = 7.44, A = 1.0 are shown 
in Fig. 3 Again we display n, Xn, Xn, d<A| 2 }, X\<p\* an d x'^ as a function 
of /i. Imthe plots we omit the data at N s = 4 since for these the transition 
point in /j, is shifted out of the /i-range of our plots due to finite volume 
effects. For the larger volumes the behavior of the observables is the same 
as for the r\ = 9.0 runs, but now the critical value [i c that marks the end of 
the Silver Blaze region is lower at /i c = 0.170(1) in lattice units. Again we 
observe scaling of the maxima of the third derivatives with the volume and 
conclude that also at rj = 7.44, A = 1.0 we find a second order transition at 
fj, c = 0.170(1). We remark that the drop of \x c = 1.146(1) in lattice units 
at rj = 9.0 to /jl c = 0.170(1) at rj = 7.44 is due to a considerable decrease 
of the lattice constant a. This implies that at rj = 7.44 we work with a 
smaller lattice spacing a and thus with a smaller physical spatial extent 
aN s . Thus we expect to see more severe finite volume effects at rj = 7.44 
than at rj = 9.0. On the other hand with the finer lattice spacing we will be 
able to go to higher temperatures (see the next section). 

To further confirm the finding of a singularity at ji c we analyzed log- 
log plots of the height of the maxima of x' n an d Xu\2 versus the 3-volume 
V = Ng. For rj = 9.0 the corresponding figure is Fig. [7] with the Nt = 100 
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Figure 3: Results at T = 0.01 for rj = 7.44 and A = 1.0. In the lhs. column 
of plots we show n, Xn and xL as a function of fi (top to bottom). In the 
rhs. column we show (|(/>| 2 ), X\<j>\ 2 an d Xup- We compare the results for 
8 3 x 100, 12 3 x 100, 16 3 x 100 and 20 3 x 100. 
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data we discuss in this section labeled by T/fj, c = 0.0087, while the r\ = 7.44 
results are shown in Fig. |] with the N t = 100 data labeled by T//i e = 0.058. 
The numerical results are represented by the symbols and we fit them with 
straight lines. The numbers next to the lines give their slope determined 
from the fits. At rj = 9.0 we find that the maxima of x'n an d xLw nicely 
follow a straight line with slopes that are close to 0.5, both for x' n an d xL\2- 
For 77 = 7.44 the finite volume effects are more severe, as discussed above, 
and we only show the data for our larger volumes, i.e, for N s = 12, 16, 20 
and 24. Some of the data show an upward curving trend and if one uses 
N s = 16, 20 and 24 for the fit, one finds slopes of 0.39(4) for the maxima of 
Xn and 0.43(11) for Xup- Although for r\ = 7.44 the data are not as clean 
due to finite volume effects (and also the temperature in units of \x c is not 
as low as for 77 = 9.0) the data indicate that also at 7? = 7.44 for the lowest 
temperature the maxima of x' n an d X L\2 diverge as V — > 00. 

The analysis implies that at zero temperature the third derivatives di- 
verge in the continuum limit and the susceptibilities develop a jump, indi- 
cating a second order transition at the end of the Silver Blaze region. In the 
first derivatives of the free energy the transition is characterized by the onset 
of a non- vanishing particle number density n at fi c and a linear response of 
the expectation value (\4>\ 2 ) to changing \i for fi > fi c . 

We conclude the discussion of the T = phenomenology with an im- 
portant cross-check. The Silver Blaze regions ends at the point where the 
chemical potential becomes equal to the mass m of the excitation it couples 
to, in other words one expects [i c = m. In order to check this we performed 
a small spectroscopy calculation on 20 3 x 40 lattices using the conventional 
representation at /i = 0. The mass was extracted from the correlator of the 
zero momentum Fourier transform of (j>, i.e., the 2-point function 

C{t) = (£(t)£(0)*> with 0(t) = _L^0(£,i). (12) 

s X 

In Fig. |4j we show our results for the effective masses m e ff(t + 1/2) = 
\nC(t)/C{t + 1) of that correlator for both, 77 = 9.0 and 77 = 7.44. The 
effective masses are represented by the symbols and the horizontal lines are 
the respective values of fi c . It is obvious that the effective masses agree 
very well with the values of /ij^] An exponential fit of the correlator for 



1 We remark that here the rj = 9.0 data are more difficult to analyze, since the large 



mass m = fi c = 1.146 leads to a very rapid exponential decay of the correlator (12 1 and 
the relative statistical error quickly overlays the signal. Thus for 77 = 9.0 the effective 
mass plateau is rather short and the errors are large. 
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Figure 4: Effective masses (symbols) for rj = 9.0 (top set of data) and 
rj = 7.44 (bottom). The effective masses are compared to the values of the 
critical chemical potential fj, c (horizontal lines) and we annotate the values 
for ji c as well as for the mass mfu obtained from an exponential fit of the 
2-point function (12). 



rj = 7.44 gives a mass of mfn = 0.171(8) which is in excellent agreement 
with [i c = 0.170(1). For rj = 9.0 we find mfn = 1.11(5) which agrees with 
He = 1.146(1) within error bars. We conclude that we properly understand 
the Silver Blaze behavior at T = . 



5 Non-zero temperature 

Having studied the transition at zero temperature we now explore how the 
transition line extends in the temperature direction of the phase diagram. 
As discussed, we work at fixed values of the couplings (again using our two 
values rj = 9.0 and rj = 7.44 at A = 1.0), and drive the temperature by 
changing Nt- This has the advantage that we always work at a fixed scale 
and do not mix temperature and discretization effects, but on the other 
hand restricts the accessible temperatures to a discrete set. For our finite 
temperature studies we use the values Nt = 100, 20, 10, 9, 8, 7, 6, 5, 4, 3 and 
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2. In lattice units the temperatures T = 1/N thus range between T = 0.01 
and T = 0.5. More interesting is the range of temperatures in units of the 
critical chemical potential at T = 0, i.e., the value [i c which we determined 
in lattice units in the previous section^} For r\ = 9.0, where fj, c = 1.146(1), 
the largest temperature we can reach with Nt = 2 is T/fi c ~ 0.436. For the 
finer rj = 7.44 lattices (fi c = 0.170(1)) we can go considerably higher, up to 
T/fic ~ 2.94, although at the price of more severe finite volume effects. 

As the temperature is increased above zero, physics cannot be completely 
independent of \x below some critical value of the chemical potential. Still, 
at least for not too large temperatures, one should observe only very weak 
/i-dependence of physical observables at sufficiently small chemical potential 
(" pseudo Silver Blaze behavior" ) . In order to define the curves in the phase 
diagram that separate the pseudo Silver Blaze behavior from a condensed 
phase characterized by strong /x-dependence, we determined the position of 
the maxima of the third derivatives of the free energy as a function of ji for 
different temperatures, using our N s = 20 lattices for rj = 9.0 and N s = 24 
for r\ = 7.44. For the r\ = 9.0 data we find that the positions of the maxima 
of x'n an d Xuia always agree exactly, while for the r\ = 7.44 data, where the 
finite size effects are more severe, we find agreement only within error bars. 

In Fig. p we show the positions of the maxima of x'n an d Xuia i n the \i-T 
plane. For The rj = 9.0 data we use filled circles connected with dashed lines 
(the plot also contains points from r] = 7.44 (triangles) which we will discuss 
below) . We express both the temperature T and the chemical potential \i in 
units of the critical value fj, c = 1.146(1) determined in the zero temperature 
runs using Nt = 100. In these units the largest temperature we can reach 
with rj = 9.0 is T/fj, c = 0.436. When turning on the temperature, the phase 
boundary curves slightly to the right reaching fj,/fj, c = 1.045 at T/fj, c = 
0.436, the highest temperature accessible for r] = 9.0. This phase boundary 
separates pseudo Silver Blaze behavior from strong /U-dependence. 

In order to assess the behavior across the phase boundary, in Fig. [6] we 
show the observables n, Xn, {\4>\ 2 } an d X|</>| 2 from the rj = 9.0 runs at Nt = 20 
as a function of \x for several temperatures (again using for T and \i units of 
fj, c ). In other words we show the observables along horizontal slices through 
the phase boundary in Fig. [5j 

While the first derivatives n and (\4>\ 2 ) were almost perfectly constant 
up to fj, c = 1.146(1) for our zero temperature runs (compare Fig. [2]), in 
Fig. [6] we now observe a deviation from constant behavior which becomes 
more prominent as the temperature is increased. However, the transition at 



2 Whenever we use fi c this refers to the critical chemical potential at zero temperature. 
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Figure 5: Phase diagram determined from the positions of the maxima of 
Xn an d Xup hi the fi-T plane. We combine data from rj = 9.0 runs for 
N s = 20 (circles connected with dashed lines) and rj = 7. 44 runs at N s = 24 
(triangles). T and // are expressed in units of /x c . 



the critical value of n is still rather pronounced, and the corresponding sus- 
ceptibilities Xn and X\$\ 2 still seem to develop a discontinuity. The position 
of the discontinuity shifts towards larger values of \i in agreement with the 
curvature of the phase boundary in Fig. [5] 

In order to analyze the order of the transition also at T > we repeat the 
analysis of the volume dependence of the maxima of the third derivatives x' n 
and Xup as function of the 3- volume V. The results of this analysis (done on 

20 3 x Nt lattices) for different temperatures are shown in Fig.[7j The symbols 
represent the numerical results and we fit them to straight lines, with the 
numbers next to the lines giving the corresponding slopes. It is obvious that 
the data are well described by a linear behavior with slopes near 0.45 for 
x' n and slopes slightly above 0.5 for Xuw- The plots show that also for our 
finite temperature values the third derivatives keep developing singularities 
in the thermodynamic limit. We conclude that for the temperature range 
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Figure 6: Results for n, Xn, (|<^| 2 ) and x\<f>\ 2 f° r V = 9.0 and N s = 20 as a 
function of \x for several temperatures corresponding to Nt = 100, 10, 5, 3 
and 2. For T and \i we use units of fi c , while n, Xn, (|0| 2 ) and X\</>\ 2 are given 
in lattice units at a fixed scale (set by r\ = 9.0, A = 1.0). 



we cover with our r] = 9.0 runs, i.e., for temperatures up to T/fj, c = 0.436 
a second order transition separates the pseudo Silver Blaze region from the 
condensed phase with strong fi dependence. 

We remark that the scaling behavior of the transition is expected to be 
different for the T = and T > cases. For T = one expects that the 
transition is of mean field type, which was indeed confirmed by Aarts in 
|12j and our dual simulation perfectly reproduces this finding (as already 
discussed in Section 3). Due to dimensional reduction for T > one expects 
a behavior in the universality class of the 3-d XY model. However, it seems 
that it is difficult to explicitly extract this behavior from a lattice simulation 
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InV InV 

Figure 7: Log-log plot of the maxima of the third derivatives of the free 
energy versus the 3-volume V = iV s 3 for tj = 9.0 and N s = 4,8, 12, 16 and 
20. We compare different temperatures corresponding to Nt = 100, 10, 5, 3 
and 2. The straight lines are linear fits to the data and the numbers near 
the lines are the respective slopes. 

with the resources available to us. This is reflected in the finding by Aarts 
|12j that the results on lattices ranging from 4 4 to 10 4 are equally well 
described by mean field behavior which illustrates that the change in the 
critical behavior must be very subtle. Also an analysis in the continuum |13| 
shows that the non-linear contribution to n versus \i — /j, c that is predicted 
for finite T is very small in the range of parameters we can access. The 
finite volume analysis in Fig. [7] is not accurate enough to extract such a 
small effect. 

An interesting question is of course where the phase boundary continues 
for larger temperatures. With the data taken at rj = 9.0 we reach a maximal 
temperature of T '/ 'fj, c = 0.436 and found the second order line to persist up to 
that value. For temperatures above that number we need to use the r\ = 7.44 
configurations where [i c = 0.170(1) in lattice units and temperatures of 
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T/nc ~ 2.94 can be reached. However, as we have already discussed, the 
rj = 7.44 data suffer from more severe finite volume effects which spoil the 
accuracy of the analysis (e.g., the peaks of the third derivatives of the free 
energy become wider and it is more difficult to exactly determine the phase 
boundary and the order of the transition). Nevertheless we try to explore 
the trend for the behavior along the phase boundary towards larger T. 

The phase diagram in Fig. [5] also contains a few points from r] = 7.44 
runs on 24 3 x Nt lattices (triangles). As for rj = 9.0 we use the positions 
of the maxima of the third derivatives x'n an d Xuia to determine the phase 
boundary. However, while for r/ = 9.0 the positions of the maxima of x'n 
and X^p always match perfectly, for rj = 7.44 this is the case only for tem- 
peratures up to T/(i, c ~ 0.5. Consequently, at rj = 7.44 we use different 
symbols for the maxima of x'n (upward pointing triangles) and Xui2 (trian- 
gles pointing downward). The two rj = 7.44 points with T '/ ' \i c < 0.5 show 
matching maxima for x' n an d Xuia an d fall on the curve set by the 77 = 9.0 
data. Those points confirm (although with larger error bars) the picture of a 
second order line up to at least T j fj, c = 0.436 from the rj = 9.0 runs. For the 
larger temperatures T/fx c > 0.5 the third derivatives x'n an d Xuia peak at 
slightly different values. However, the positions of the maxima agree within 
error bars and we attribute the small discrepancy to the more severe finite 
volume effects for the rj = 7.44 ensembles. The corresponding observables 
are shown in Fig. [8] where we plot n, Xn, (|</>| 2 ) an d X\<j>\ 2 f° r V = 7.44 as a 
function of \x for several temperatures. 

In Fig. [9] we repeat the finite volume study of Fig. [7] now for rj = 7.44. 
We show the log-log plot of the maxima of x'n an d Xuw versus the spatial 
volume V = N^. It is obvious that here the data do not as nicely fall on a 
straight line as in the corresponding plots for 77 = 9.0 shown in Fig. [7j We fit 
only the data points for the three largest volumes and as in Fig. [7] denote the 
slopes next to the lines. It is quite interesting to observe that the slope for the 
maxima of Xup remains at values close to 0.5 up to the largest temperatures, 
indicating that the susceptibility X|</>| 2 manages to develop a discontinuity 
in the thermodynamic limit also for our highest temperatures. The stronger 
finite volume effects at 77 = 7.44 make the analysis more difficult for x'n- 
The slopes for x n are smaller for the smaller volumes but show a tendency 
towards pointing upwards for larger V, indicating second order behavior for 
the //-derivatives also up to the highest temperatures we can access. 

One expects that the second order line persists for large T due to the 
breaking of the U(l) symmetry at the transition to the condensed phase. In 
Fig. 10 we attempt a plausibility test for this scenario: We plot the height of 
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Figure 8: Results for n, Xn, (M 2 ) and x\<j>\ 2 f° r V = 7.44 and N s = 24 as a 
function of \x for several temperatures corresponding to Nt = 100, 20, 10, 8 
and 6. For T and n we use units of fj, c , while n, Xn, (\</>\ 2 ) and X|</>| 2 are given 
in lattice units at a fixed scale (set by rj = 7.44, A = 1.0). 



the maxima of the third derivatives x' n an d Xuia (both made dimensionless 
with suitable powers of (i c ) as a function of the inverse temperature in units 
of ji c . We find that the maxima are non-zero at finite T and only for 1/T — > 
the curves reach 0. Furthermore the maxima for the r\ = 9.0 runs are 
larger than those at rj = 7.44 as expected from the larger volume at rj = 9.0. 
Although this analysis is only an exploratory test at a fixed volume, it 
supports the existence of a phase transition throughout the phase diagram. 

We conclude our presentation of the finite temperature results with the 
discussion of observables along two vertical cuts in the phase diagram Fig. [5] 
Using the rj = 9.0 data at N s = 20, in Fig. fTTlwe show the observables n, 
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Figure 9: Log-log plot of the maxima of the third derivatives of the free 
energy versus the 3-volume V = Nf for rj = 7.44 and N s = 12, 16, 20 and 24. 
We compare different temperatures corresponding to Nt = 100, 10, 6 and 3. 
The straight lines are linear fits to the data (using only N s = 16, 20, 24) and 
the numbers near the lines are the respective slopes. 



Xm (H 2 ) and x\</)\ 2 as a function of T/fi c and compare two different values 
of the chemical potential. The vertical cut at /u//u c = 1.012 hits the phase 
boundary at T//j, c ~ 0.2, while the cut at fi/fi c = 0.977 continues all the 
way to T / n c = without entering the condensed phase (compare Fig. [5]). 
This behavior is clearly reflected in the observables: For the n/n c = 0.977 
data the particle number density n and the susceptibility Xn both drop to 
zero as the temperature is decreased to T j \x c = 0, as expected for a value of 
the chemical potential in the Silver Blaze domain, i.e., fi < /i c . The vacuum 
expectation value X\4>\ 2 drops towards 0.237, while the corresponding suscep- 
tibility X\<t>\ 2 vanishes as T j \x c — >■ 0. For the cut at fj,/ fj, c = 1.012 the behavior 
is quite different: When decreasing the temperature down to T/fi c ~ 0.2, 
where the cut hits the phase boundary of the condensed phase, the sus- 
ceptibilities show a drastic increase an then remain essentially constant for 
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Figure 10: Size of the maxima of the third derivatives x'n an d xUp (both 
made dimensionless with suitable powers of fi c ) as a function of the inverse 
temperature in units of \i c . We compare the results for rj = 9.0 on 20 3 x Nt 
lattices and r] = 7.44 on 24 3 x N t . 



T/[A C < 0.2. The particle number density n has a non-vanishing value at 
T/nc = and (\4>\ 2 ) settles at a value of 0.285. The plots nicely demonstrate 
the temperature driven condensation of the system for fi/fi c = 1.012 and its 
absence at /x/// c = 0.977. 



6 Summary and discussion 

In this work we have studied a charged scalar <ji 4 field at finite density using 
the lattice formulation. The complex action problem of the standard form 
was overcome by mapping the partition sum to a flux representation where 
only real and non-negative terms appear. The system was updated with a 
generalized worm algorithm and we studied thermodynamical observables 
as a function of temperature and chemical potential. 

For zero temperature the system undergoes a second order phase transi- 
tion at some critical value \i c that marks the end of the Silver Blaze region. 
For fj, < fi c the observables are independent of /z, while above a strong /in- 
dependence sets in. The critical value [i c could very cleanly be identified 
from the narrow peaks of two 3rd derivatives of the free energy and in a 
small spectroscopy calculation we showed that fj, c coincides with the mass of 
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the scalar field. A finite volume analysis established the second order nature 
of the transition. 

Increasing the temperature we found that the critical line bends towards 
larger values of fi. For temperatures up to T j \x c ~ 0.5 we could reliably es- 
tablish the continuation of second order behavior for finite T. At small but 
non-zero temperatures we found only a very weak dependence of observables 
on the chemical potential for values of fi below the critical value (pseudo Sil- 
ver Blaze behavior) . At the critical line we saw an abrupt change of physical 
observables marking the transition into the condensed phase characterized 
by strong //-dependence. When studying our observables as a function of the 
temperature we could nicely illustrate the temperature driven condensation 
for /j,/ He = 1.012 and the absence of condensation at (jl/ (jl c = 0.977. 
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For the analysis of temperatures T/fj, c > 0.5 we had to rely on finer 
lattices (rj = 7.44) which, however, show considerably more severe finite size 
effects. For small temperatures the results match the results from our coarser 
(77 = 9.0) lattices very well, but for T//j, c > 0.5 some of the observables - 
the ones related to the particle number density - were seriously affected 
by finite volume effects and thus gave inconclusive results about the nature 
of the transition between pseudo Silver Blaze behavior and the condensed 
phase with strong /i-dependence. The observables related to derivatives of 
the free energy with respect to the mass parameter, i.e., (|^| 2 ), X\4>\ 2 an d 
Xu,i2> on the other hand indicate that a second order line could persist to 
temperatures considerably higher than T/fx c ~ 0.5. Lattice studies on large 
and fine lattices will be necessary to settle this issue. 
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Appendix: Pseudo-code for the worm algorithm 



In this appendix we provide pseudo-code for the initialization of the worm 
and its subsequent steps. The worm algorithm that updates the variables 
k XjV £ Z is set up in a fixed background of the variables l XjU generated 
with local Metropolis sweeps (see Section [3]) . The starting point of a worm 
will be denoted by xo, the position of its head by x, and A will be used 
for the set of all lattice sites. A is an equally distributed random variable 
with A £ {— 1, +1} that determines whether the worm tries to increase or 
decrease the flux along its path. By v £ {±1, ±2, ±3, ±4} we denote the 
direction of the next attempted step of the worm, and in case v is negative 
[y = we use the convention k x> _\ v \ = ^ X _Q \ u \ f° r the- assignment of 

the corresponding link variables. As in the main text we use the abbrevi- 
ation f x = Y,y[\ k xA + \K-V,v\ + 2(i z ,i/ + lx-v,u)], and by f x we denote the 
same expression but evaluated with the variable k X]V from the trial offer. Fi- 
nally, rand ( ) denotes a random generator for uniformly distributed random 
numbers in the interval [0, 1). 

Starting the worm 

while worm has not started 

choose A £ { — 1,+1} randomly 

choose xo £ A randomly 

choose v £ {±1, ±2, ±3, ±4} randomly 

kxo,p — k X q V + A 

if 1^0,^1 ^ |^xo,y| then 

p = exp (sign(u) fiA5\ u \ A ) x (\k XOjl/ \ +l Xo ,u)~ X Wif^^y 1 

else 

p = exp (sign(v)(iA6\ v \ >4 ) x (\k XQ:U \ +l Xo ,„) x Wif^+^y 1 

endif 



if rand() < p then 

X = Xo + V 

worm has started 

endif 



end while 
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Running the worm 
while x / xo 

choose v £ {±1, ±2, ±3, ±4} randomly 

x = x + P 

if |fc XjI/ | > |fc XjI/ | then 

p = exp (s^n(^) / uA(5| i ,| i 4)x(|fc :Eii ,|+/ a: ^) _ xW(f x )xW(f x+ z 

else 

p = exp (sign(u)/j,A6^ A )x(\k XOtU \+l XOtV )xW(f x )xW(f x+i ;) 
endif 

if x = x then pi — pxW(f Xo ) 
if rand() <p then 

X i X 

endif 
end while 
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